

function [objw,wresL,gammaw_hat]=MMw(q, p, y, x, w, toler, gamma)
  
 % Returns the QR estimator for a quantile regression model under weighted
 % residuals
 % The QR estimator is computed using a MM (Majorize-Minimize) algorithm
 % from Hunter Lange (2000) JCGS
 
 %******************* Input Variables **************************************

 % q : quantile level
 % p : dimension of the covariate (including intercept)
 % y : dependent variable
 % x : covariate (must include intercept)
 % w : weights used in the random weighting bootstrap method
 % toler : tolerance parameter for the precision of the estimation
 % gamma : initial set of value for the intercept and slope coefficients
 
 %******************* Output Variables **************************************
 
 % objw : optimized objective function
 % wresL : residuals of the estimation
 % gammaw_hat : vector storing the estimates of the intercept and slope
 % coefficients
 
  n=length(y);
  
  %******** MM tolerance and precision variables ****************************
  
  tn=toler/n;
  e0=-tn/log(tn);
  epsilon=(e0-tn)/(1+log(e0)); %  epsilon is the smoothing parameter.  
  % As suggested in Hunter Lange paper, it
  %  is chosen to roughly satisfy -epsilon(log epsilon)=toler/n using a
  %  Newton-Raphson step above.
  
  %******** Algorithm initialisation *******************
  
  iteration=0;
  change=realmax;
  res=y-x*gamma; % Compute residuals for current gamma.
  wres=w.*res; % weight the residuals
  weights=1./(epsilon+abs(wres));
  v=1-2*q-wres.*weights;
  lastsurr=-sum(wres.*v)-sum((1-2*q).*wres); % Last value of surrogate (up
                                             % to a constant).
  
  %******** Algorithm loop *******************       
  
  while change > toler 
	iteration = iteration + 1;
	J=x.*w(:,ones(p,1));  % J is the first differential matrix.

       	step = -inv(J'*(weights(:,ones(p,1)).*J)) * (J' * v);
	    gamma=gamma+step;
       	res=y-x*gamma;
        wres=w.*res;
        surr = sum(wres.^2.*weights) + sum((4*q - 2).*wres);
        
        %  Now do the step-halving procedure to ensure a decrease in the
        %  value of the surrogate function.
        while surr > lastsurr
       		step=step/2;
       		gamma=gamma-step;
       		res=y-x*gamma;
            wres=w.*res;
       		surr = sum(wres.^2.*weights) + sum((4*q - 2).*wres);
       	end

       	change=lastsurr-surr;
	    weights = 1./(epsilon+abs(wres));
	    v = 1-2*q - wres.*weights;
	    lastsurr = -sum(wres.*v)-sum((1-2*q).*wres);
  end
  gammaw_hat=gamma;
  wresL=wres;
  objw = sum(q.*wres)-sum(wres(wres<0));

end
